library(tidyverse)
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(RColorBrewer)
library(treemapify)
source("Custom_function.R")
Read Bidoup - Nui Ba data with read.xlsx() function
#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# Bidoup - Nui Ba plot
BDplot <- read.xlsx(filename, sheet = "BD_census2")
We practice with 2ha data, in which BD-1 (1ha) is concentrated (100x100m) and BD-2 (1ha) is taken from 25 scattered subplots.
Use n_distinct() counts the number of distinct value in
column: family, genus, scientificName. And
group_by() to count at each plot.
BDplot |>
group_by(PlotID) |>
summarise("family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName),
"stem" = n())
## # A tibble: 2 × 5
## PlotID family genus species stem
## <chr> <int> <int> <int> <int>
## 1 BD-1 46 70 124 12207
## 2 BD-2 53 94 169 11575
We usually group trees into 3 size-class of DBH:
small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm
BDplot <- BDplot |>
mutate(
dbh_cm = dbh / 10,
# dbh value in mm
gDBH = case_when(dbh_cm >= 20 ~ ">=20",
dbh_cm >= 5 ~ "[5-20)",
TRUE ~ "[1-5)")
)
head(BDplot)
We count the number of stems at each diameter class and each plot
(stemByDBH <- BDplot |>
group_by(PlotID, gDBH) |>
summarise(n = n()) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n)))
## `summarise()` has grouped output by 'PlotID'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 4
## # Groups: PlotID [2]
## PlotID gDBH n freq
## <chr> <chr> <int> <dbl>
## 1 BD-1 [1-5) 9736 0.798
## 2 BD-2 [1-5) 9478 0.819
## 3 BD-1 [5-20) 2057 0.169
## 4 BD-2 [5-20) 1749 0.151
## 5 BD-1 >=20 414 0.034
## 6 BD-2 >=20 348 0.03
Figure
ggplot(BDplot, aes(factor(gDBH), fill = factor(PlotID))) +
geom_bar(position = "dodge") +
theme_bw(14) +
labs(x = "DBH (cm)", y = "Nb of stem", fill = "PlotID")
We can also count the number of species at each diameter class.
(
SpByDBH <- BDplot |>
group_by(PlotID, gDBH) |>
summarise(
"family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName)
)
)
## `summarise()` has grouped output by 'PlotID'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 5
## # Groups: PlotID [2]
## PlotID gDBH family genus species
## <chr> <chr> <int> <int> <int>
## 1 BD-1 >=20 28 35 53
## 2 BD-1 [1-5) 45 69 120
## 3 BD-1 [5-20) 42 61 96
## 4 BD-2 >=20 33 45 64
## 5 BD-2 [1-5) 53 91 160
## 6 BD-2 [5-20) 48 74 121
Figure
# Convert to longer table
SpByDBH <- SpByDBH |>
tidyr::pivot_longer(-c(PlotID,gDBH), names_to = "rank", values_to = "n")
Plot BD-1
SpByDBH |> filter(PlotID=="BD-1") |>
ggplot(aes(x = gDBH, y = n, fill = rank)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = n, label = n),
position = position_dodge(width = 0.9),
vjust = -0.5) +
scale_fill_viridis_d() +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Rank")
Plot BD-2
SpByDBH |> filter(PlotID=="BD-2") |>
ggplot(aes(x = gDBH, y = n, fill = rank)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = n, label = n),
position = position_dodge(width = 0.9),
vjust = -0.5) +
scale_fill_viridis_d() +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Rank")
BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
Draw a blank plot with draw_blank_plot_20ha()
function.
p <- draw_blank_plot_20ha()
p
Add add tree to map
p + geom_point(data = BDplot.BD1,
aes(x = gx, y = gy, size = dbh_cm, color = gDBH),
alpha = 0.7) +
xlim(80, 180) + ylim(400, 500) +
labs(color = "DBH class", size = "DBH(cm)")
For big tree
BDplotA <- BDplot.BD1 |> filter(gDBH == ">=20")
p + geom_point(data = BDplotA,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For medium tree
BDplotB <- BDplot.BD1 |> filter(gDBH == "[5-20)")
p + geom_point(data = BDplotB,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For small tree
BDplotC <- BDplot.BD1 |> filter(gDBH == "[1-5)")
p + geom_point(data = BDplotC,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For a particular species, px: Pinus krempfii
## So do phan bo toan bo cay trong o mau
BDplotSP <- BDplot.BD1 |> filter(scientificName == "Pinus krempfii")
p + geom_point(data = BDplotSP,
aes(x = gx, y = gy, size = dbh_cm, color = gDBH),
alpha = 0.7) +
xlim(80, 180) + ylim(400, 500) +
labs(color = "DBH class", size = "DBH(cm)")
We can also try with a detail DBH class
BDplot.BD1 <- BDplot.BD1 |>
mutate(
dbh_cm = dbh / 10 , # dbh in mm
gDBH2 = case_when(
dbh_cm >= 100 ~ ">=100",
dbh_cm >= 90 ~ "[90-100)",
dbh_cm >= 80 ~ "[80-90)",
dbh_cm >= 70 ~ "[70-80)",
dbh_cm >= 60 ~ "[60-70)",
dbh_cm >= 50 ~ "[50-60)",
dbh_cm >= 40 ~ "[40-50)",
dbh_cm >= 30 ~ "[30-40)",
dbh_cm >= 20 ~ "[20-30)",
dbh_cm >= 10 ~ "[10-20)",
TRUE ~ "[1-10)")
)
stemByDBH2 <- BDplot.BD1 |>
count(gDBH2) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n))
stemByDBH2
## gDBH2 n freq
## 1 [1-10) 11169 0.915
## 2 [10-20) 624 0.051
## 3 [20-30) 246 0.020
## 4 [30-40) 123 0.010
## 5 [40-50) 35 0.003
## 6 [70-80) 3 0.000
## 7 >=100 2 0.000
## 8 [50-60) 2 0.000
## 9 [60-70) 1 0.000
## 10 [80-90) 1 0.000
## 11 [90-100) 1 0.000
stemByDBH2 |>
ggplot(aes(x = gDBH2, y = n)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5) +
#scale_fill_manual(values = c("grey","white"))+
theme_bw(14) +
labs(x = "DBHClass (cm)", y = "Number of stem")
The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:
relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)
\[IVI={\left(relative.density + relative.basal.area \right)}\]
Basal area (BA) of each species can be calculated from diameter measured: \(BA = pi * R^2\)
Relative basal area (RBA): \(RBA = BA / sum(BA)\)
Relative density (RD) = number of individuals of each species / total
BDplot <- BDplot |>
mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)
BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
BDplot.BD2 <- BDplot |> filter(PlotID == "BD-2")
IVI.1 <- calcIVI(BDplot.BD1)
head(IVI.1)
## # A tibble: 6 × 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus pierrei 966 6.74 12.7 7.91 20.6 1
## 2 Polyosma nhatrangensis 1175 2.33 4.40 9.63 14.0 2
## 3 Syzygium rubicundum 726 3.21 6.06 5.95 12.0 3
## 4 Pinus krempfii 89 5.28 9.96 0.729 10.7 4
## 5 Castanopsis echinocarpa 539 3.05 5.75 4.42 10.2 5
## 6 Lasianthus saprosmoides 1076 0.236 0.444 8.81 9.26 6
IVI.2 <- calcIVI(BDplot.BD2)
head(IVI.2)
## # A tibble: 6 × 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus coinhensis 357 3.94 7.84 3.08 10.9 1
## 2 Castanopsis echinocarpa 483 3.27 6.51 4.17 10.7 2
## 3 Polyosma nhatrangensis 705 1.31 2.60 6.09 8.69 3
## 4 notdet 883 0.127 0.253 7.63 7.88 4
## 5 Rhaphiolepis poilanei 436 1.07 2.13 3.77 5.90 5
## 6 Lithocarpus pierrei 262 1.40 2.79 2.26 5.05 6
Figure
figureIVI(IVI.1)
figureIVI(IVI.2)
biodivIndex(IVI.1)
## Richness Simpson Shannon Evenness
## 1 124 0.9598452 3.765398 0.7811574
biodivIndex(IVI.2)
## Richness Simpson Shannon Evenness
## 1 169 0.9756952 4.200136 0.8187561
Above-ground biomass (AGB)
Below-ground biomass (BGB)
Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.
AGB=Total aboveground biomass (kg)
D=DBH (cm)
H=Height (m)
\(ρ\)=Wood density (g/cm3)
BA=Basal area (cm2)
| Equation | |
| Tropical forest DRY | \(AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)\) |
| \(AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\) | |
| Tropical forest MOIST | \(AGB = exp(-2.977 + ln(ρ * D^2 * H))\) |
| \(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\) |
library(BIOMASS)
## For more information on using BIOMASS, visit https://umr-amap.github.io/BIOMASS
## Using temporary cache
## It is recommended to use a permanent cache to avoid to re-download files on each session.
## See function createCache() or BIOMASS.cache option.
woodensity <- getWoodDensity(
genus = BDplot.BD1$genus,
species = BDplot.BD1$scientificName,
family = BDplot.BD1$family
)
## The reference dataset contains 16467 wood density values
## Your taxonomic table contains 124 taxa
head(woodensity)
BDplot.BD1$meanWD <- woodensity$meanWD
\(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)
BDplot.BD1 <- BDplot.BD1 %>%
mutate(dbh_cm = dbh/10,
AGB = meanWD * exp(-1.499 + 2.148*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))
(AGB <- sum(BDplot.BD1$AGB)/1000) #Mg/ha
## [1] 469.1828
The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level
(BGB = 0.275 * AGB)
## [1] 129.0253
Total biomass of 1 hecta
(TB <- AGB + BGB)
## [1] 598.208
The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.
#Carbon stock (ton/ha)
(TC <- TB * 0.47)
## [1] 281.1578
# CO2 (ton)
(CO2 <- TC * 44/22)
## [1] 562.3155
1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.
In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)
1 hecta forest at Bidoup - Nui Ba national park ~ 2800 $
(CO2 * 5)
## [1] 2811.578
LITHCA : Lithocarpus eucalyptifolius
CASTEC: Castanopsis echinocarpa
SYMPBA: Symplocos banaensis